function phiGrad = cfdComputeGradientGaussLinear0(phi)
%--------------------------------------------------------------------------
%
%  Written by the CFD Group @ AUB, Fall 2018
%  Contact us at: cfd@aub.edu.lb
%==========================================================================
% Routine Description:
%   This function calculates gradient based on green gauss method
%   The values are interpolated to the faces using a weighted factor with no
%   correction for non-conjunctionality, hence the 0 digit in the function 
%   name. It takes as input the phi element array and returns the element 
%   gradient array.
%--------------------------------------------------------------------------

% Check field type (by size)
theSize = size(phi,2);
if theSize==3
    theNumberOfComponents = 3;
else
    theNumberOfComponents = 1;
end

% Get mesh
theMesh = cfdGetMesh;

% INTERIOR FACES contribution to gradient
iFaces = 1:theMesh.numberOfInteriorFaces;
% get the list of owners and neighbours for all interior faces
owners_f = cfdGetOwnersSubArrayForInteriorFaces;
neighbours_f = cfdGetNeighboursSubArrayForInteriorFaces;

% get the surface vector of all interior faces 面法向量
Sf = cfdGetFaceSfSubArrayForInterior;
% get the geometric factor for all interior faces （用于从element向face上插值）
faceWeights = cfdGetFaceWeightsForInterior;

% Initialize phiGrad Array
phiGrad = zeros(theMesh.numberOfElements+theMesh.numberOfBElements,3,theNumberOfComponents);

% Accumulate face fields and areas
for iComponent=1:theNumberOfComponents    
    % compute the linear interpolation of phi into all interior faces
    phi_f = faceWeights.*phi(neighbours_f,iComponent) + (1-faceWeights).*phi(owners_f,iComponent);
    % loop over faces and add contribution of phi face flux to the owner and 
    % neighbour elements of the faces
    for iFace=iFaces
        % 注意对于owners是正值，对于neighbors是负值，因为面的法向量Sf值的符号是owners正neighbors负
        % 面的法向量Sf的大小为这一个面的面积，这里算梯度是用了Gauss公式将 grad(phi)
        % 的体积分化成
        phiGrad(owners_f(iFace),:,iComponent)     = phiGrad(owners_f(iFace),:,iComponent)     + phi_f(iFace)*Sf(iFace,:);
        phiGrad(neighbours_f(iFace),:,iComponent) = phiGrad(neighbours_f(iFace),:,iComponent) - phi_f(iFace)*Sf(iFace,:);
    end    
end

% BOUNDARY FACES contribution to gradient
iBElements = theMesh.numberOfElements+1:theMesh.numberOfElements+theMesh.numberOfBFaces;
% get the list of elements owning a boundary face
owners_b = cfdGetOwnersSubArrayForBoundaryPatchFaces;
% get the boundary values of phi
phi_b = phi(iBElements,:);
% get the surface vector of all boundary faces
Sb = cfdGetFaceSfSubArrayForBoundaryPatchFaces;
% loop over all boundary faces and add phi flux contribution to gradient
for iComponent=1:theNumberOfComponents
    for k=1:theMesh.numberOfBFaces
        phiGrad(owners_b(k),:,iComponent) = phiGrad(owners_b(k),:,iComponent) + phi_b(k)*Sb(k,:);
    end
end

% Get Average Gradient by dividing with element volume
% get volume of all elements in mesh
volumes = cfdGetVolumesForElements;
% compute the average gradient by dividing the sum of all phi fluxes for an 
% element by the volume of the element
for iComponent=1:theNumberOfComponents
    for iElement =1:theMesh.numberOfElements
        phiGrad(iElement,:,iComponent) = phiGrad(iElement,:,iComponent)/volumes(iElement);
    end
end

% Set cfdBoundary Gradient equal to associated element Gradient
% 再次注意，element field除了element外，还存储了与boundary faces一一对应的
% boundary element（它俩实际就是一个东西，值相等，只不过是将face值存在了element的数组里）
phiGrad(iBElements,:,:) = phiGrad(owners_b,:,:);

end